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Abstract 



In the present work fournontrivial stages of electrokinetic instability are identified by di- 
rect numerical simulation (DNS) of the full Nernst-Planck-Poisson-Stokes (NPPS) sys- 
tem: i) The stage of the influence of the initial conditions (milliseconds); ii) ID self-similar 
O ■ evolution (milliseconds-seconds); iii) The primary instability of the self-similar solution 
(seconds); iv) The nonlinear stage with secondary instabilities. The self-similar charac- 
i-£h 1 ter of evolution at intermediately large times is confirmed. Rubinstein and Zaltzman 
instability and noise-driven nonlinear evolution to over-limiting regimes in ion-exchange 

CN ! membranes are numerically simulated and compared with theoretical and experimental 

> ■ 

VjD ■ predictions. The primary instability which happens during this stage is found to arrest 
OO ' 

^ 1 self-similar growth of the diffusion layer and specifies its characteristic length as was first 
experimentally predicted by Yossifon and Chang (PRL 101, 254501 (2008)). A novel 
q ! principle for the characteristic wave number selection from the broadbanded initial noise 

is established. 

Introduction. Pattern formation in dissipative systems has a long and intriguing his- 

X- 

tory in a number of disciplines. It is generally associated with nonlinear effects which 
lead to qualitatively new phenomena that do not occur in linear systems [1]. Problems 
of electro-kinetics has recently attracted a great deal of attention due to the rapid de- 
velopment of micro-, nano- and biotechnologies. A curious and fascinating electrokinetic 
instability and pattern formation in an electrolyte solution between semi-selective ion- 
exchange membranes under a DC potential drop was theoretically predicted in [2,3] and 
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experimentally confirmed in [4], including direct experimental proof [5-7]. Both theory 
and experiments show that these instabilities are reminiscent of Rayleigh-Benard con- 
vection and Benard-Marangoni thermo-convection, but are much more complicated from 
both the physical and mathematical points of view. 

In the present work for the first time nontrivial stages of the noise-driven electroki- 
netic instability are obtained by direct numerical simulation (DNS) of the full Nernst- 
Planck-Poisson-Stokes (NPPS) system. Small-amplitude initial room disturbances have 
the following stages of evolution: i) A transient stage of the influence of the ini- 
tial conditions, t = 0(X 2 D /AD) (milliseconds); ii) A ID self-similar stage of evolution, 
X 2 D /4D <C t <C L 2 /4Z) (milliseconds-seconds); iii) A stage in which instability becomes 
manifest and arrests any self-similar growth and selects a diffusion layer characteris- 
tic length, t = 0(L 2 /AD) (seconds); iv) A nonlinear stage with secondary instabilities, 
t > L 2 /AD. 

Formulation. A binary electrolyte between perfect semi-selective ion-exchange mem- 
branes is considered. Notations with tilde are used for the dimensional variables, as op- 
posed to their dimensionless counterparts without tilde. The diffusivities of cations and 
anions are assumed to be equal, D + = D~ = D. The following characteristic values are 
taken to make the system dimensionless: distance between the membranes, L; dynamic 
viscosity, jl; thermodynamic potential, $o = RT/F; bulk ion concentration at initial 
time, cq. Then, the characteristic time and characteristic velocity are scaled, respec- 
tively, by L 2 /D and D/L. Here F is Faraday's constant; R is the universal gas constant; 
T is Kelvin temperature; jl is the dynamic viscosity of the fluid; e is the permittivity of 



the medium; and the dimensional Debye length is taken as = y eRT / F 2 cq. 

Electro-convection is then described by the transport equations for the concentrations 
for positive and negative ions, c + , c~ , the Poisson equation for the electric potential, 
and the Stokes equation for a creeping flow, presented in dimensionless form: 



dc* 

z/ 2 V 2 $ = c" - c+ VP - V 2 U = kV 2 $V$, V • U = 



— + U ■ Vc* = ±V ■ (^V^) + V 2 ^, ^ 



with the boundary conditions at the membrane surfaces, y = and 1, a fixed concen- 
tration of positive ions, no flux condition for negative ions, and a fixed drop of potential 



and velocity adhesion: 

<9$ dc~ 

c + = p, - c ~— + — = 0, $ = 0fory = 0, $ = AVfory=l, U = 0, (2) 
ay ay 

with the characteristic electric current j at the membrane surface y = 0, j = c + d§/dy + 
dc + /dy. 

The problem is described by three dimensionless parameters: the drop of potential, 
AV; the dimensionless Debye length, v = Xd/L, and the coupling coefficient between 
hydrodynamics and electrostatics, k = e^/jlD (the last is fixed for a given electrolyte 
solution). The dependence on the wall concentration, p, for over-limiting regimes is 
practically absent [3] and, hence, p is not included in the mentioned parameters. 

The typical bulk concentration of the aqueous electrolytes varies in the range Co = 
1 -z- 10 3 mol/m 3 ; the potential drop is about AV = -r 5 V; the absolute temperature 
can be taken to be T = 300° K; the diffusivity is about D = 2 • 10 -9 m 2 /s; the distance 
between the electrodes L is of the order of 0.5 -r 1.5 mm; the concentration p on the 
membrane surface must be much larger than c and it is usually taken within the range 
from 5co to 10co- The dimensional Debye layer thickness Ad is varied within the range 
from 0.5 to 15 nm depending on the concentration c - 

The DNS of the system (1), without any simplification, is implemented by apply- 
ing the Galerkin pseudo-spectral r-method. A periodic domain along the membrane 
surface allows the utilization of a Fourier series, exp(nkx), in the x-direction. Cheby- 
shev polynomials, T m (y), are applied in the transverse direction, y. The accumula- 
tion of zeros of these polynomials near the walls allows of properly resolving the thin 
space charge regions. Eventually, the physical variables are presented in the form 
^ = Em En ^nm(^ni(!/) exp (nb). A special method is developed to integrate the 
system in time. The details of the numerical scheme will be presented elsewhere [8]. The 
basic wave number k is connected with the membrane length / as k = 2n/l. Then, k = 1 
and I = 2ir are taken such that the dimensional length of the membrane / corresponds to 
the experimental one in [5]. The problem is solved for v = 10~ 3 and 5- 10 -4 , dimensionless 
drop of potential AV — -r 200 (AV^ = -r 5 V), and for a typical value of k, k — 0.2; 
p = 5. The results of the full-scale numerical investigation will be accompanied by a 

3 



discussion of the main stages of the evolution along with a comparison with simplified 
solutions [9-11]. 

Self-similar stage of evolution. Adding initial conditions makes the statement (l)-(2) 
complete. This system is a strongly nonlinear one and a change in the initial conditions 
can change the attractor (see, for example, famous experimental work by Gollub and 
Swinney [12] where several attractors, depending on the initial conditions, were found 
for flow between two rotating cylinders). Hence what must be taken are conditions that 
are natural from the experimental view-point. From this view-point, when the drop of 
potential between the membranes is zero, there are only two thin double-ion layers near 
the membrane surfaces and the ion distribution is uniform outside these layers. This 
specifies the initial conditions. Then, a small-amplitude white-noise spectrum should be 
superimposed on the bulk ion concentrations, c + = c~ = 1, 

t = : c + — 1 = white noise(x), c~ — 1 = white noise(x), (3) 

non-uniform along the membrane. Because of the smallness of the amplitude the behavior 
at initial times can be assumed one-dimensional (ID), as if the right-hand side of (3) were 
zero, and, as a consequence of this one- dimensionality, any convection and non-uniformity 
along x can be discarded from consideration at initial times: U = 0, d/dx = 0. 

When the potential drop is turned on, an extended space charge (ESC) and diffusion 
regions are developed near the bottom membrane. There is an interval of intermediate 
times, \ 2 D /AD < t < L 2 /4L>, or in dimensionless form, jis 2 <^t<^. \, when the problem 
does not have a static characteristic size: the double-ion layer is too small and the distance 
between the membranes is too large. On dimensional grounds, the only diffusion length 
which can give a proper characteristic size that includes time and is, hence, dynamical, 
is - and the behavior of the solution has to be self-similar. These intermediate 

times vary from milliseconds to seconds (for details see [9-11]). 

The self-similar rescaling of the variables, along with the ID assumption, turns the 
PDE system (l)-(2) into the ODE system with the new dimensionless spatial variable 
rj = yj \/ 4Dt and a small parameter e — Ad/ V 4Dt. A new dimensionless Debye length 
e, electric current J and drop of potential A$ in self-similar variables are connected with 
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Figure 1: (Left) Distribution of the charge density p in space for several moments of time, 
AV = 100 and v = \d/L — 0.001; (a) — numerical solution of (1)~(3), (b) — self-similar 
solution. The numerical points for v = 0.0005, and several values of AV shrink into one 
self-similar VC curve (shown by the solid line). 



the corresponding old variables, v, j and AV, by the relations 



e(t) 



J(t) = 2j(t)Vt, A$(t) = AV 



(4) 



The two-parameter family of self-similar solutions for the parameters A<J> and e is 
found in [9-11] and now is compared with the DNS of (l)-(3). In Fig. 1 (left), such a 
comparison for the charge distribution p{rj) = c + (rj) — c~{rf) is presented for different mo- 
ments of time where the parameters of the self-similar solutions are recalculated according 
to (4). In Fig. 1 (right) universal voltage-current characteristics J versus AF = eA$ are 
shown by a solid line while triangles, stars and circles are gathered for several fixed drops 
of potential between the membranes, AV, for intermediate times. The numerical points 
for all AV shrink into the universal self-similar VC curve. A rather good correspondence 
justifies a self-similar behavior at intermediate times. 

Instability of the self-similar solution and the length scale selection. The initial stage 
of evolution is one-dimensional and self-similar, but room disturbances (3) for the super- 
critical case, AV > AV, will grow and eventually manifest themselves, destroying the 
self-similarity and making the solution two-dimensional. This specifies the next step 
of evolution. When the amplitudes of the harmonics are small, their behavior can be 
considered independently. Imposing on the ID self-similar solution a perturbation of the 




Figure 2: Marginal stability curve at AV = 40 and v = 0.001. The time-dependent coor- 
dinates are inverse Debye length 1/e and wave number a. Each straight line corresponds 
to a perturbation with a given wave number k\ a point on the line is moving with time 
towards the stability curve and eventually crosses it. 
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Figure 3: Wave number selection for the room disturbances; darker regions correspond 
to larger amplitudes. Broadband initial noise is transformed into a sharp band near 
k ~ 32 4- 34, dashed line corresponds to the prediction of stability theory, fc* = 35, 
AV = 50 and u = 0.001. 
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Figure 4: Typical evolution of the average electric current (j)/ju m (solid line) and current 
amplitude j max — j m ; n (dashed line) (left) and diffusion layer thickness (right) from room 
disturbances for v = 0.001 and AV = 60. Three time intervals are exhibited; point 1 
specifies the instability selected length scale. 

form R(r), t) exp (iax) and linearizing (l)-(3) with respect to R turns this system into 
a linear PDE system with time-dependent coefficients. Restricting ourselves to neutral 
perturbations, d/dt = 0, transforms this PDE system into an ODE system in t]. Solution 
of the eigenvalue problem for this system gives the marginal stability conditions [9]. 

The parametric dependence of the Debye length and wave number for the self-similar 
solution in time completely changes the interpretation of the marginal stability curves [11]. 
It is instructive to present the stability curve at a fixed drop of potential AV in coordinates 
of the inverse Debye length l/e = 2y/i/ v versus the wave number a = kvje = 2k\/i, both 
referred to the self-similar length and, so, time-dependent, see Fig. 2. Each straight 

line in the figure corresponds to a monochromatic wave with a fixed wave number, k = aL; 
the larger is k, the steeper is the line. Inside the marginal stability curve a perturbation 
is growing, but outside — decaying. A point on the straight line corresponding to some 
initial perturbation is moving with time towards the marginal stability curve just because 
of the parametrical time dependence of axes. At the same time, the amplitude of the 
disturbance at the beginning is always decaying but after crossing the marginal curve at 
some point, it starts to grow. The nose critical point with coordinates (l/£*,o:*) first 
reaches the marginal stability curve and, hence, first loses its stability. Triangles in the 
figure are taken from our DNS for monochromatic perturbations. There is a rather good 
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quantitative agreement between the numerical simulations and the stability theory for 
self-similar solutions. 

Room disturbances at the linear stage of evolution can be considered as a superposition 
of individual harmonics. Their evolution always starts in the stable region and filters a 
broadband initial noise (3) into a sharp band of wave numbers. It stands to reason that we 
may assume that the wave number selection is determined by the critical points (1/e*, «*) 
with the characteristic wave number fc* = a^e^jv. This principle is different from the 
well-known selecting principle of maximum growing mode and is verified by comparison of 
the results of the DNS with the BC's (3): the comparison for typical conditions presented 
in Fig. 3 shows the validity of this novel principle. 

The relation A(t) = 2/ (j(t)) between the diffusion layer thickness A and the electric 
current averaged along the membrane length /, (J) = ~ j(x,t)dx, is accepted. For the 
self-similar solution the diffusion layer thickness should be proportional to the square 
root of time, A(t) ~ y/i, and, respectively, (j(t)} ~ l/\J\t). The typical evolution of 
the electric current (J) and of the diffusion layer thickness A, presented in Fig. 4, exhibit 
three characteristic time intervals. In the interval II, for times of intermediate magnitude, 
A is proportional to \ft and the behavior is one- dimensional and self-similar, while for 
shorter times, in the time interval I, as well as for longer times in the interval III, this 
self-similarity is violated. In region I, the influence of the initial conditions is important, 
while in region III, the growth of A is arrested by instability. The current amplitude 
Jmax — Jmin characterizes instability and is shown by the dashed line in Fig. 4 (left); it 
becomes significant at the point 1 which separates regions II and III and specifies the 
characteristic length of the diffusion layer. Further calculations show that this length 
selected by instability is not changed much by subsequent nonlinear processes. 

In Fig. 5, a qualitative comparison with the experiments of Yossifon and Chang [7] 
is presented (for details see [10]). The time dependence of the thickness of the diffusion 
layer A is taken in dimensional variables. Excluding very short times of establishment, 
up to 8.9 s of evolution, the solution in zone I follows the self-similar law proportionally 
to Vt (dash-and-dot line); then the self-similar process is violated by the loss of stability 
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Figure 5: Time dependence of the thickness of the diffusion layer A in dimensional 
variables for AV = 2.5 V, L = 1 mm, D = 5 • 10~ 9 m 2 /s and v = 0.001 [10]. (1) Self- 
similar solution, (2) numerical solution, I is the self- similarity range, and II is the stability 
loss by the trivial solution and the appearance of electroconvection. The corresponding 
experimental dependence [7] is in the inset. 

in zone II. According to the numerical experiment, the transition to the electroconvection 
mode takes place in the range 6 -v- 10 s. This corresponds to the value of 8 s observed in 
physical experiments [7]. 

Nonlinear stages of evolution and secondary instability. The next stages of evolu- 
tion are nonlinear processes with eventual saturation of the disturbance amplitude. For 
a small super-criticality, AV" — AV*, the nonlinear saturation leads to steady periodic 
electro-convective vortices along the membrane surface; with increasing super-criticality 
the attractor can be described as a structure of periodically oscillating vortices; with 
further increase in super-criticality, the behavior eventually becomes chaotic in time and 
space. 

The evolution of the normalized to jn m average electric current and current amplitude 
of perturbation are presented in Fig. 6. For small super-criticality this evolution results 
in steady vortices, Fig. 6 (left), and for larger super-criticality — in chaotic behavior, 
Fig. 6 (right). 
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Figure 6: Evolution of the normalized average electric current (j)/jn m (solid line) and 
current amplitude (j max — jmm)/iiim (dashed line) for v = 0.001 k = 0.2 for small super- 
criticality, AV = 30 (left) and for large super-criticality, AV = 50 (right). 




Figure 7: Snapshots of the charge density p(x, y, t) at different times. Darker regions 
correspond to larger charge densities (for parameters, see previous figure). 
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Snapshots of the charge density p(x,y,t) for the above mentioned conditions at dif- 
ferent times are shown in Fig. 7. Darker regions correspond to larger densities of space 
charge. There is a rather sharp boundary between the ESC region and the diffusion 
region. Point 1 corresponds to the initial small-amplitude broadbanded white noise. A 
linear stability mechanism filters it into practically periodic disturbances with a wave 
number fc*, point 2. The boundary between the ESC region and the diffusion region be- 
haves sinusoidally: its minimum corresponds to the maximum of the charge density and 
vice versa. With increasing time, there is a fusion of neighboring waves with correspond- 
ing reducing of their number. The sinusoidal profile changes to the spike-like one where 
regions of small charge in the spikes are jointed by thin flat regions of large space charge, 
moment 5, left. Besides fusion of the spikes there is also their creation so that eventually 
some equilibrium of their number is established, moment 5, right. Movie visualization for 
charge density evolution can be found in [13] (AV = 50, v = 10 -3 and k = 0.2, chaotic 
behavior) . 

Let us consider the physical mechanisms of the secondary instability. The sharp 
boundary between the ESC and diffusion regions has small charge densities near its 
humps and large charge densities in neighboring regions around its hollows, either for 
periodic or spike-like coherent structures. These regions with larger charge are trying to 
expand (charge of the same sign is repelling) and this can lead to the disappearance of 
some of the humps — with eventual coarsening. Such instability is illustrated in Fig. 8 
(left). There is also another physical mechanism which results in the opposite effect. 
If a flat region between two humps is long enough, it suffers from primary Rubinstein 
and Zaltsman electroconvective instability, see Fig. 8 (right). Points 1, 2 and 3 are 
nucleation points of future humps. They grow and eventually form new humps. This last 
mechanism is applicable only for spike-like structures with flat neighboring regions, but 
not for sinusoidal waves. 

Conclusions. The results of the DNS reported in this paper refer to realistic exper- 
imental parameters that may be realized in lab experiments. The following stages of 
evolution are identifed: i) a stage in which the initial conditions have considerable in- 
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Figure 8: Two physical mechanisms of the secondary instability, (a) Coalescence of two 
humps with low charge density with following coarsening, (b) Instability of flat region 
between humps with following creation of new humps. 

fluence (milliseconds); ii) ID self-similar evolution (milliseconds-seconds); iii) primary 
instability of the self-similar solution (seconds); iv) a nonlinear stage with secondary in- 
stabilities. This work can be viewed as a new step in the understanding of the instabilities 
and pattern formation in micro- and nano-scales. 

The authors are grateful to I. Rubinstein, M. Z. Bazant, and H.- C. Chang for fruitful 
discussions. This work was supported in part by Russian Foundation for Basic Research 
(project no. 11-08-00480-a). 
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Figure Captions 

Figure 1. (Left) Distribution of the charge density p in space for several moments of time, 
AV = 100 and v = Xd/L = 0.001; (a) — numerical solution of (l)-(3), (b) - 
self-similar solution. The numerical points for v = 0.0005, and several values of 
AV shrink into one self-similar VC curve (shown by the solid line). 

Figure 2. Marginal stability curve at AV = 40 and v = 0.001. The time-dependent co- 
ordinates are inverse Debye length 1/e and wave number a. Each straight line 
corresponds to a perturbation with a given wave number k\ a point on the line is 
moving with time towards the stability curve and eventually crosses it. 

Figure 3. Wave number selection for the room disturbances; darker regions correspond to 
larger amplitudes. Broadband initial noise is transformed into a sharp band near 
k pa 32 4- 34, dashed line corresponds to the prediction of stability theory, k* = 35, 
AV = 50 and v = 0.001. 

Figure 4. Typical evolution of the average electric current {j)/jn m (solid line) and current 
amplitude j max — j min (dashed line) (left) and diffusion layer thickness (right) from 
room disturbances for v = 0.001 and AV = 60. Three time intervals are exhibited; 
point 1 specifies the instability selected length scale. 

Figure 5. Time dependence of the thickness of the diffusion layer A in dimensional variables 
for AV = 2.5 V, L = 1 mm, D = 5 • 10~ 9 m 2 /s and v = 0.001 [10]. (1) Self- 
similar solution, (2) numerical solution, I is the self-similarity range, and II is the 
stability loss by the trivial solution and the appearance of electroconvection. The 
corresponding experimental dependence [7] is in the inset. 

Figure 6. Evolution of the normalized average electric current (j) /jn m (solid line) and current 
amplitude (j max — Jnun)/jiim (dashed line) for v = 0.001 k = 0.2 for small super- 
criticality, = 30 (left) and for large super-criticality, AV = 50 (right). 

Figure 7. Snapshots of the charge density p(x, y, t) at different times. Darker regions corre- 
spond to larger charge densities (for parameters, see previous figure). 
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Figure 8. Two physical mechanisms of the secondary instability, (a) Coalescence of two humps 
with low charge density with following coarsening, (b) Instability of flat region 
between humps with following creation of new humps. 
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